#NOTE: based on Ward & Ahlquist 2018, chapter 3
rm(list=ls()); gc()
library(foreign); library(MASS); library(ggplot2)

inv.logit <- function(x) {
  exp(x)/(1+exp(x))
}

set.seed(1234)

setwd("")
getwd()

# load data 
cohesion <- read.dta("datacross.dta")

# Variables
cohesion$dyadid <- as.integer(cohesion$dyadid)
cohesion$MKactive <- as.factor(cohesion$MKactive)
cohesion$frag_postMK1 <- as.numeric(cohesion$frag_postMK1)
cohesion$irregular <- as.numeric(cohesion$irregular2)
cohesion$rebsupnum <- as.numeric(cohesion$rebsupnum)
cohesion$excluded <- as.numeric(cohesion$excluded)
cohesion$prevactive<- as.numeric(cohesion$prevactive)
cohesion$fightcapnum<- as.numeric(cohesion$fightcapnum)
cohesion$fightcontnum<- as.numeric(cohesion$terrcontnum)
cohesion$intensityfirstyear <- as.numeric(cohesion$intensityfirstyear)
cohesion$dyadid <- as.integer(cohesion$dyadid)
cohesion$exclpop <- as.numeric(cohesion$exclpopn)
cohesion$e_v2x_neopat <- as.numeric(cohesion$e_v2x_neopat)
cohesion$dyadsfirstyearconflict <- as.numeric(cohesion$dyadsfirstyear)

# fitting the model:
withcov <- glm(frag_postMK1 ~ MKactive + frag_preMK1 + irregular2 + rebsupnum + excluded +
                 prevactive + fightcapnum + terrcontnum + e_v2x_neopat + dyadsfirstyear + exclpop + 
                 intensityfirstyear, family=binomial(link="logit"), data=cohesion, na.action=na.omit)

summary(withcov)


x.lo <- c(1, #intercept
          0, # MKactive == 0
          median(cohesion$frag_preMK1, na.rm=TRUE),
          median(cohesion$irregular2, na.rm=TRUE),
          median(cohesion$rebsupnum, na.rm=TRUE),
          median(cohesion$excluded, na.rm=TRUE),
          median(cohesion$prevactive, na.rm=TRUE),
          median(cohesion$fightcapnum, na.rm=TRUE),
          median(cohesion$terrcontnum, na.rm=TRUE),
          median(cohesion$e_v2x_neopat, na.rm=TRUE),
          median(cohesion$dyadsfirstyear, na.rm=TRUE),
          median(cohesion$exclpop, na.rm=TRUE),
          median(cohesion$intensityfirstyear, na.rm=TRUE))
X.lo <- matrix(x.lo, nrow=length(x.lo), 1)

x.hi <- c(1, #intercept
          1, # MKactive == 1
          median(cohesion$frag_preMK1, na.rm=TRUE),
          median(cohesion$irregular2, na.rm=TRUE),
          median(cohesion$rebsupnum, na.rm=TRUE),
          median(cohesion$excluded, na.rm=TRUE),
          median(cohesion$prevactive, na.rm=TRUE),
          median(cohesion$fightcapnum, na.rm=TRUE),
          median(cohesion$terrcontnum, na.rm=TRUE),
          median(cohesion$e_v2x_neopat, na.rm=TRUE),
          median(cohesion$dyadsfirstyear, na.rm=TRUE),
          median(cohesion$exclpop, na.rm=TRUE),
          median(cohesion$intensityfirstyear))

X.hi <- matrix(x.hi, nrow=length(x.hi), 1)

B.tilde <- mvrnorm(1000, coef(withcov), vcov(withcov)) #1000 draws of coefficient vectors
s.lo <- inv.logit(B.tilde %*% X.lo) #matrix of predicted probabilities
s.hi <- inv.logit(B.tilde %*% X.hi) #matrix of predicted probabilities

## median
s.lo<-apply(s.lo, 2, quantile, c(0.025, 0.5, .975), na.rm = FALSE)
s.hi<-apply(s.hi, 2, quantile, c(0.025, 0.5, .975), na.rm = FALSE)

s.hi
s.lo

# matrix of results:
ResMatrix <- matrix(NA, 2,3)
colnames(ResMatrix) <- c("lower","median","upper")
rownames(ResMatrix) <- c("NoTargeting", "YesTargeting")
ResMatrix[1,1] <- s.lo[1,]
ResMatrix[1,2] <- s.lo[2,]
ResMatrix[1,3] <- s.lo[3,]
ResMatrix[2,1] <- s.hi[1,]
ResMatrix[2,2] <- s.hi[2,]
ResMatrix[2,3] <- s.hi[3,]
ResMatrix

# figure 1
pdf("figure1.pdf", width=5.5, height=5.5) 
ggplot(as.data.frame(ResMatrix), aes(x = c(0,1))) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.5, width=.01) +
  geom_point(aes(y = median), size = 1.5) + theme_bw() + theme(legend.text=element_text(size=12)) +
  labs(x = "State-led Collective Targeting", y = "Probability(Fragmentation)") +
  ylim(-0.05, 0.55) + scale_x_continuous(breaks = c(0,1))
dev.off()
